Maximum-likelihood estimation prevents unphysical Mueller matrices 
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We show that the method of maximum-hkelihood estimation, recently introduced in the context of 
quantum process tomography, can be apphed to the determination of MueUer matrices characterizing 
the polarization properties of classical optical systems. Contrary to linear reconstruction algorithms, 
the proposed method yields physically acceptable Mueller matrices even in presence of uncontrolled 
experimental errors. We illustrate the method on the case of an unphysical measured Mueller matrix 
taken from the literature. 

OCIS codes: 000.3860, 260.5430, 290.0290. 



In the mathematical description of both polarized light 
and two-level quantum systems (or quhits, in the lan- 
guage of quantum information) , there are many analogies 
and common tools. For example, the Poincare sphere 
for classical polarization and the Bloch sphere for two- 
level quantum systems |^ are, in fact, the same mathe- 
matical object. Although the classical concepts and tools 
were introduced well before the quantum ones, the lat- 
ter were developed independently of the former. Thus, 
many well established results in classical polarization op- 
tics have been "rediscovered" in the context of quantum 
optics and quantum information 3] . Interestingly, the in- 
verse process of borrowing results from quantum to clas- 
sical optics has started only recently |j, Q . 

In this Letter we give a contribution to this inverse 
process by pointing out a connection between quantum 
process tomography (QPT) Q and classical polarization 
tomography (CPT). Specifically, we show that the re- 
cently introduced maximum-likelihood (ML) method for 
the estimation of quantum processes [ilSjiSjj can be suc- 
cessfully applied to the determination of classical Mueller 
matrices. In the conventional approach to CPT, Mueller 
matrices are estimated from the measurement data by 
means of linear algorithms 10]. However, such recon- 
structed Mueller matrices often fail to be physically ac- 
ceptable We show that this problem is avoided by 
using the maximum-likelihood method which allows to 
include in a natural manner the "physical-acceptability" 
constraint. Thus, thanks to a "quantum tool", we solve 
an important issue that has been long debated in the 
classical literature 0, ll^l Q| . This is in particular im- 
portant in view of the present interest in CPT, e.g., for 
medical and astronomical imaging. 

To begin with, we give first a qualitative description of 
the connection between QPT and CPT. At the heart of 
this connection lies the well known mathematical equiv- 
alence (isomorphism) between the density matrix p de- 
scribinga two-level quantum system and the coherency 
matrix [l| J describing the classical polarization state of 
a light beam [HIi^: 



J/Tr J. 



(1) 



is p. A quantum process that transform an input state 
p'" in an output state can be described by a linear 
superoperator Q : — Qp^^. Analogously, a classical 
linear optical process (as, e.g., an elastic scattering pro- 
cess), can be described by a 4 x 4 matrix Ai such that 
jout _ or, in explicit representation. 



k,l 



M^.mJ'Cu i,j,fc,/e{0,l}. (2) 



In the same way as the reconstruction of Q is the goal 
of QPT, the estimation of the elements AAij.ki from the 
measurement data is the goal of CPT. However, in the 
common practice, instead of the complex matrix ele- 
ments Adij^ki one wants to determine the 16 real elements 
M^jy, {fj,,iy = 0, . . . , 3) of the associated Mueller matrix 
M. In the ML approach, the estimated elements of M 
are found to be the most likely to yield the measured 
data. In what follows we show how to find them. 

In a classical polarization tomography experiment, the 
measurement data are collected following the scheme 
shown in Fig. 1. An input light beam is prepared in a 
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J is an Hermitian, positive semidefinite 2x2 matrix, as 



FIG. 1: Classical polarization tomography (CPT) scheme. 



pure polarization state represented by the coherency ma- 
trix J'j^y and sent through the optical system S where 
it is transformed in the output beam represented by 
"^(b) ~ ^"^(b)- '^^^ estimation strategy is to retrieve 
information on the system from a series of polarization 
tests on the output states J°^^ obtained from distinct in- 
put states J"^y {b — 1, . . . , B). A polarization filter V{a) 
that allows the passage of light with specific polariza- 
tion labelled by the index a, [a — 1, . . . , A), provides for 
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the polarization tests. Finally, the intensity lab of the 
beam after the filter is recorded. If we prepare B differ- 
ent input states and we perform A polarization tests per 
each output state, then such a CPT experiment will have 
Ax B outcomes corresponding to all measured intensities 
{lab}, ia^l,...,A,b^l,...,B). 

However, in a ML approach, which is a probabilistic 
method, one deals with relative rather than absolute in- 
tensities. Therefore, from the data set {lab} we must 
extract the relative intensities (or measurement frequen- 
cies) fab = lab/h, where lb is the intensity of the in- 
put light beam. Since by definition < fab < 1 and 
fab + fa'b £ 1, whcrc tt, a' label two mutually orthogonal 
polarization tests, then fab provides an experimental es- 
timation of the theoretical probability Pab for obtaining 
a nonzero output intensity with polarization a when the 
input beam is prepared in the polarization state b. The 
theoretical probabilities {pab} can be written in terms of 
the input and output states as [3| 



Pab(M) 



Tr[n(,)j; 



out] 
(a)-'(b) \ 



Tr j; 



(3) 
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where we denoted with H^a) — e(a) (g) e|^j the 2x2 pro- 
jection matrix representing the action of the polarization 
filter V{a) oriented along the (possibly complex) unit vec- 
tor e(a). Since in a CPT experiment the input beam is 
always prepared in a pure polarization state, it can be 
represented by a projection matrix n(fc-) as J^'^Jj = /bn(6), 

where h — Tr J^'"^ is the intensity of the beam. Then we 
can rewrite 



PabiM) = Tr[n(„)A^n(fc)] 



(4) 



where from now on we assume, without loss of generality, 
/{, = !. At this point, having measured the frequencies 
{fab} and having calculated the probabilities {pab}, the 
sought matrix A4 can be obtained by maximizing a like- 
lihood function -C[A^] defined as 



C[M] = J2f-bHPab{M)] 



a,b 

^4blnTr[n(„)X%)] 

a,b 



(5) 



where Pab{-M) > for any physically acceptable process. 

Equation © is the first main result of this pa- 
per. It contains both experimental (lab) and theoreti- 
cal [pab{M)] quantities. Now, we demonstrate that it 
is possible to impose the condition Pafc(A^) > before 
the maximization operation, in such a way that the esti- 
mated Mueller matrix is automatically physically accept- 
able. After a lengthy but straightforward calculation, it 
is possible to show that the matrix M can be written in 
terms of an Hermitian matrix H as llSl 



M 



0,3 



(m) ^ ^Hl ' 



(6) 



where {((p,)} are the elements of the standard basis in 
C^^^. By substituting Eq. © into Eq. © we obtain 
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Pab{M) = ^i?^^Tr[n(a)e(^)n(f,)e(^) 



(7) 



where the superscript T indicates the transposed matrix. 
The probabilities Pab{M.) as written in Eq. Q can still 
be negative, because only H matrices associated with 
physically acceptable Mueller matrices can guarantee the 
condition pab > 0. However, we know from the Mueller 
matrix theory that the H matrix associated with a physi- 
cally acceptable Mueller matrix must be positive semidef- 
inite [l^ . It is well known that any positive semidefinite 
matrix can be written in terms of its Cholesky decompo- 
sition as 



H 



^Tr(CCt)' 
where C is a lower triangular matrix 



C 



(8) 
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(9) 



composed by 16 real parameters hk, {k = 1, . . . , 16), and 
we fixed the normalization of H by setting Mqo — 1. 
Then, after substituting Eq. © into Eq. {Tj), the maxi- 
mum of C can be found by using a standard maximiza- 
tion algorithm 19]. The search for the maximum is per- 
formed in the real 16-dimensional space of parameters 
{/ife}. Once the optimal set of values {h°^^ , . . . , h°Q^ } 
that maximize C has been found, this can be used in Eq. 
(O to obtain the corresponding i/°P*. Finally, the ele- 
ments of the sought physically acceptable Mueller matrix 
can be computed as 



M^. = Tr{F°P*[a(^)®ay} 



(10) 



where {o'(^)} are the normalized Pauli matrices jlS] . This 
is our second main result. A Mueller matrix M deter- 
mined in this way represents the answer to the question: 
which physically acceptable Mueller matrix is most likely 
to yield the measured data? " 

The rest of the paper is devoted to the illustration of 
the theory outlined above, by applying it to a realistic 
case. We have chosen from the current literature the 
following Mueller matrix which was already shown 
to be physically unacceptable: 



M' = 



( 0.7599 0.0295 0.1185 -0.0623\ 

0.0384 0.5394 0.0282 -0.1714 

0.124 -0.012 0.6608 0.2168 

\-0.0573 -0.1811 -0.1863 0.4687 / 



(11) 
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From M' we calculated the (normalized) associated Her- 
mitian matrix H' which is not positive semidefinite since 
it has one negative eigenvalue: 

diagi?' = {1.6671, 0.2950, 0.2330, -0.1951}. (12) 

By using Eq. Q, we generated a set of 36 "fake- 
measured" data as 



0,3 



(13) 



where we selected both the input beam (represented by 
!!({,)) and the polarization filter (represented by '^{^a)) 
from the set of 6 pure polarization states labelled as hor- 
izontal (H) and vertical (V), oblique at 45° and oblique 
at 135°, right (RHC) and left (LHC). The so obtained 
36 numbers represent our "experimental" data set. Ob- 
viously, from these numbers one could generate back M' 
with the conventional linear algorithm. It may worth to 
note that the three pairs of polarization states we have 
chosen, are the ones usually employed in CPT_[1|, and 
correspond to three mutually unbiased basis [2l| often 
utihzed in QPT. We used the MATHEMATICA 5.1 func- 
tion FindMaximum to maximize C Tho run this func- 
tion it is necessary to furnish a set of initial values for 
the parameters {hi, . . . , hi^} to be estimated. We found 
convenient to proceed in the following way: we made 
first a polar decomposition of H' to obtain the positive 

semidefinite matrix H" = 2^/ H'H''' /TtV H'H'\ then 
we obtain the initial values from the Cholesky decompo- 
sition of H" . Finally, after maximization, we obtained 
the maximum-likelihood estimation of as 



/ 0.7599 
0.0372 
0.1208 
\-0.0554 



0.0257 
0.5285 
-0.0001 
-0.0572 



0.1206 
0.0001 
0.6184 
-0.1794 



-0.0576\ 
-0.0496 
0.1920 
0.4822 / 



(14) 

As expected, this matrix is indeed a physically acceptable 
Mueller matrix, as the eigenvalues of its associated 
matrix are all non-negative: 

diagiJ^'L = {1.6344, 0.2341, 0.1315, 0}. (15) 

A visual inspection show that M and differs only 

by a little amount. This was expected since we choose 
an initial Mueller matrix M' that is not very unphysical 
(only one negative eigenvalue ). A quantitative estima- 
tion of the difference between M' and can be given 
by calculating their relative Frobcnius distance \u\ 



\\M -M 



ML I 



\\M + M^^\\ 



0.072, 



(16) 



about 7%. This confirms the quality of our approach 
even with sparse data set (only 36 values). 



In conclusion, we have shown that it is possible to ap- 
ply the maximum-likelihood method, initially developed 
for quantum process tomography, to the classical prob- 
lem of Mueller matrix reconstruction. Moreover, we have 
shown that this method has the benefit to produce always 
physically acceptable Mueller matrices as the most likely 
matrices which yield the measured data. 
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which indicates that the average relative difference be- 
tween corresponding matrix elements of M and is 
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